Full-time away home signal

Lag analysis at a state-level

Use most correlated shift time to build linear models

\[y = \beta_{0} + \beta_{1}(\sum I(A_i = 1)) + \beta_{2}x_1 + \beta_{3}x_2\]

# select best lag number based on intial exploration
case.shifted.days.spearman <- 37 # based spearman correlation
doc.visit.shifted.days.spearman <- 62 # based spearman correlation for doctor visit
cum.death.shifted.days.spearman <- 55
death.shifted.days.spearman <- 55

case.vec <- c("confirmed_7dav_incidence_num", "confirmed_7dav_cumulative", "confirmed_7dav_cumulative_prop")
doc.vec <- c("smoothed_cli", "smoothed_adj_cli")
cum.death.vec <- c("deaths_7dav_cumulative_num")
death.vec <- c("deaths_7dav_incidence_num")

# Make two copies for shifting the data
selected.df <- intervention_mobility_case
factored_data <- intervention_mobility_case
# Change the data by shifting the covariates
# Shift the case count column vector by the specified shift time
factored_data <- shiftDays(selected.df, factored_data, case.shifted.days.spearman, case.vec)

# Shift the  doctor column vector by the specified shift time
factored_data <-shiftDays(selected.df, factored_data, doc.visit.shifted.days.spearman, doc.vec)
  
# Shift the death count column vector by the specified shift time
factored_data <-shiftDays(selected.df, factored_data, cum.death.shifted.days.spearman, cum.death.vec)

# Shift the cum death count column vector by the specified shift time
factored_data <-shiftDays(selected.df, factored_data, death.shifted.days.spearman, death.vec)
  
# We specifically look at emergency declaration
factored_data%>% 
  mutate(EmergDec.duration = cumsum(EmergDec)) %$% 
  lm(full_time_work_prop ~ EmergDec.duration + smoothed_cli+smoothed_adj_cli ) %>% 
  summary()
## 
## Call:
## lm(formula = full_time_work_prop ~ EmergDec.duration + smoothed_cli + 
##     smoothed_adj_cli)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0169158 -0.0058953  0.0006305  0.0053490  0.0235429 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.515e-02  1.367e-03  25.719  < 2e-16 ***
## EmergDec.duration  7.255e-05  1.670e-05   4.345 2.16e-05 ***
## smoothed_cli      -4.762e-03  8.335e-04  -5.714 3.72e-08 ***
## smoothed_adj_cli   5.119e-03  1.035e-03   4.945 1.55e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.007749 on 212 degrees of freedom
##   (458 observations deleted due to missingness)
## Multiple R-squared:  0.278,  Adjusted R-squared:  0.2678 
## F-statistic: 27.21 on 3 and 212 DF,  p-value: 6.297e-15
# We specifically look at emergency declaration
factored_data%>% 
  mutate(EmergDec.duration = cumsum(EmergDec)) %$% 
  lm(full_time_work_prop ~ EmergDec.duration + smoothed_cli+smoothed_adj_cli+confirmed_7dav_cumulative_prop+deaths_7dav_cumulative_num) %>% 
  summary()
## 
## Call:
## lm(formula = full_time_work_prop ~ EmergDec.duration + smoothed_cli + 
##     smoothed_adj_cli + confirmed_7dav_cumulative_prop + deaths_7dav_cumulative_num)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.016885 -0.005988  0.000543  0.005137  0.023104 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     3.848e-02  4.761e-03   8.083 6.17e-14 ***
## EmergDec.duration               2.387e-05  7.231e-05   0.330    0.742    
## smoothed_cli                   -4.780e-03  8.525e-04  -5.607 6.85e-08 ***
## smoothed_adj_cli                5.492e-03  1.221e-03   4.500 1.16e-05 ***
## confirmed_7dav_cumulative_prop -2.570e-06  6.187e-06  -0.415    0.678    
## deaths_7dav_cumulative_num      9.238e-07  1.555e-06   0.594    0.553    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.007872 on 198 degrees of freedom
##   (470 observations deleted due to missingness)
## Multiple R-squared:  0.2535, Adjusted R-squared:  0.2346 
## F-statistic: 13.45 on 5 and 198 DF,  p-value: 2.712e-11

Model with time only (weekends included)

Model with shifting the covariates (weekends included)

Model without shifting the covariates (weekends included)

intervention.lm <- intervention_mobility_case %>% 
  mutate(EmergDec.duration = cumsum(EmergDec)) 

lm.fit.no.lag <- lm(full_time_work_prop ~ EmergDec.duration + smoothed_cli+smoothed_adj_cli, data =intervention.lm
) 

intervention.lm$predlm <- c(rep(NA, nrow(intervention.lm) - length(predict(lm.fit.no.lag))), predict(lm.fit.no.lag))

intervention.lm%>% 
  mutate(policy.duration = cumsum(EmergDec), EmergDeclaration = as.factor(EmergDec)) %>% 
  ggplot(aes(x = time_value, y = full_time_work_prop, color = EmergDeclaration)) +
  geom_point() + 
  geom_line(aes(x = time_value, y = predlm, colour="fitted value"), size = 1)+
  labs(title = "Covariates selected WITHOUT most correlated number of shift")

Weekend effects

We suspect that the mobility signal is lower than usual during the weekend.

intervention_mobility_case$weekday <- weekdays(as.Date(intervention_mobility_case$time_value)) 

p <- ggplot(intervention_mobility_case, aes(x=weekday, y=full_time_work_prop)) + 
  geom_boxplot()
p

Re-plot the regression line after dropping weekends

Regressing on time only

Covariates shifted

## 
## Call:
## lm(formula = full_time_work_prop ~ EmergDec.duration + smoothed_cli + 
##     smoothed_adj_cli)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0208421 -0.0026275 -0.0001302  0.0029902  0.0131662 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.976e-02  1.022e-03  38.918  < 2e-16 ***
## EmergDec.duration  5.185e-05  1.244e-05   4.169 5.16e-05 ***
## smoothed_cli       7.097e-04  7.307e-04   0.971    0.333    
## smoothed_adj_cli  -2.629e-04  8.398e-04  -0.313    0.755    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.004876 on 150 degrees of freedom
##   (328 observations deleted due to missingness)
## Multiple R-squared:  0.4393, Adjusted R-squared:  0.4281 
## F-statistic: 39.18 on 3 and 150 DF,  p-value: < 2.2e-16

Covariates not shifted

## 
## Call:
## lm(formula = full_time_work_prop ~ policy.duration + smoothed_cli + 
##     smoothed_adj_cli)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0305809 -0.0058044 -0.0005214  0.0055933  0.0242504 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       5.576e-02  1.248e-03  44.671  < 2e-16 ***
## policy.duration   5.029e-05  1.346e-05   3.736 0.000245 ***
## smoothed_cli      1.837e-03  9.248e-04   1.986 0.048441 *  
## smoothed_adj_cli -4.955e-03  1.208e-03  -4.104 5.99e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.008831 on 194 degrees of freedom
##   (284 observations deleted due to missingness)
## Multiple R-squared:  0.1633, Adjusted R-squared:  0.1504 
## F-statistic: 12.62 on 3 and 194 DF,  p-value: 1.427e-07

What will happen if we add more covidcast signals as covariates?

Adding covidcast signals

## 
## Call:
## lm(formula = full_time_work_prop ~ EmergDec.duration + smoothed_cli + 
##     smoothed_adj_cli + confirmed_7dav_incidence_num + confirmed_7dav_cumulative + 
##     confirmed_7dav_cumulative_prop + deaths_7dav_incidence_num + 
##     deaths_7dav_cumulative_num)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0197354 -0.0031287 -0.0000458  0.0030361  0.0123833 
## 
## Coefficients: (1 not defined because of singularities)
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     4.811e-02  4.510e-03  10.667  < 2e-16 ***
## EmergDec.duration              -9.391e-05  7.241e-05  -1.297  0.19682    
## smoothed_cli                    7.370e-04  7.269e-04   1.014  0.31239    
## smoothed_adj_cli               -1.310e-03  9.767e-04  -1.342  0.18189    
## confirmed_7dav_incidence_num    1.026e-06  3.302e-07   3.106  0.00230 ** 
## confirmed_7dav_cumulative      -1.019e-08  1.179e-08  -0.864  0.38886    
## confirmed_7dav_cumulative_prop         NA         NA      NA       NA    
## deaths_7dav_incidence_num       1.037e-04  3.222e-05   3.218  0.00161 ** 
## deaths_7dav_cumulative_num      1.677e-06  1.294e-06   1.296  0.19700    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.004806 on 138 degrees of freedom
##   (336 observations deleted due to missingness)
## Multiple R-squared:  0.4433, Adjusted R-squared:  0.4151 
## F-statistic:  15.7 on 7 and 138 DF,  p-value: 4.839e-15

Adding other policy duration

# Try to add other intervention covariates

factored_data.without.weekend %$%  
  lm(full_time_work_prop ~EmergDec.duration +StayAtHomeDuration+PublicMaskDuration+SchoolCloseDuration+GathRestrictDuration+BarRestrictDuration+NEBusinessCloseDuration+ RestaurantRestrictDuration+  smoothed_cli+smoothed_adj_cli+confirmed_7dav_incidence_num+confirmed_7dav_cumulative+confirmed_7dav_cumulative_prop+deaths_7dav_incidence_num+deaths_7dav_cumulative_num+ SchoolCloseDuration) %>% 
  summary()
## 
## Call:
## lm(formula = full_time_work_prop ~ EmergDec.duration + StayAtHomeDuration + 
##     PublicMaskDuration + SchoolCloseDuration + GathRestrictDuration + 
##     BarRestrictDuration + NEBusinessCloseDuration + RestaurantRestrictDuration + 
##     smoothed_cli + smoothed_adj_cli + confirmed_7dav_incidence_num + 
##     confirmed_7dav_cumulative + confirmed_7dav_cumulative_prop + 
##     deaths_7dav_incidence_num + deaths_7dav_cumulative_num + 
##     SchoolCloseDuration)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0197354 -0.0031287 -0.0000458  0.0030361  0.0123833 
## 
## Coefficients: (8 not defined because of singularities)
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     4.811e-02  4.510e-03  10.667  < 2e-16 ***
## EmergDec.duration              -9.391e-05  7.241e-05  -1.297  0.19682    
## StayAtHomeDuration                     NA         NA      NA       NA    
## PublicMaskDuration                     NA         NA      NA       NA    
## SchoolCloseDuration                    NA         NA      NA       NA    
## GathRestrictDuration                   NA         NA      NA       NA    
## BarRestrictDuration                    NA         NA      NA       NA    
## NEBusinessCloseDuration                NA         NA      NA       NA    
## RestaurantRestrictDuration             NA         NA      NA       NA    
## smoothed_cli                    7.370e-04  7.269e-04   1.014  0.31239    
## smoothed_adj_cli               -1.310e-03  9.767e-04  -1.342  0.18189    
## confirmed_7dav_incidence_num    1.026e-06  3.302e-07   3.106  0.00230 ** 
## confirmed_7dav_cumulative      -1.019e-08  1.179e-08  -0.864  0.38886    
## confirmed_7dav_cumulative_prop         NA         NA      NA       NA    
## deaths_7dav_incidence_num       1.037e-04  3.222e-05   3.218  0.00161 ** 
## deaths_7dav_cumulative_num      1.677e-06  1.294e-06   1.296  0.19700    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.004806 on 138 degrees of freedom
##   (336 observations deleted due to missingness)
## Multiple R-squared:  0.4433, Adjusted R-squared:  0.4151 
## F-statistic:  15.7 on 7 and 138 DF,  p-value: 4.839e-15
# Predict the mobility
new.pred <- factored_data.without.weekend %$%  
  lm(full_time_work_prop ~EmergDec.duration +StayAtHomeDuration+PublicMaskDuration+SchoolCloseDuration+GathRestrictDuration+BarRestrictDuration+NEBusinessCloseDuration+ RestaurantRestrictDuration+  smoothed_cli+smoothed_adj_cli+confirmed_7dav_incidence_num+confirmed_7dav_cumulative+confirmed_7dav_cumulative_prop+deaths_7dav_incidence_num+deaths_7dav_cumulative_num+ SchoolCloseDuration)%>%
  predict()

# Pad the fitted values with NA 
factored_data.without.weekend$predlm <- c(rep(NA, nrow(factored_data.without.weekend) - length(new.pred)), new.pred)

# Plot the graph
factored_data.without.weekend %>%
  ggplot(aes(x = time_value, y = full_time_work_prop, color = EmergDeclaration)) +
  geom_point() + 
  geom_line(aes(y = predlm, colour="fitted value"), size = 1) +
   labs(title = "All covariates selected WITH most correlated number of shift (weekends dropped)")